Elliptic and triangular flow in event-by-event (3+l)D viscous hydrodynamics 
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We present results for the elliptic and triangular flow coefficients vi and V3 in Au+Au collisions 
at y/s = 200AGeV using event-by-event 3+fD viscous hydrodynamic simulations. We study the 
effect of initial state fluctuations and finite viscosities on the flow coefficients V2 and V3 as functions 
of transverse momentum and pseudo-rapidity. Fluctuations are essential to reproduce the measured 
centrality dependence of elliptic flow. We argue that simultaneous measurements of «2 and vs can 
determine n/s more precisely. 
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Fluctuating initial conditions for hydrodynamic sim- 
ulations of heavy-ion collisions have been argued to be 
very important for the exact determination of collective 
flow observables and to describe specific features of multi- 
particle correlation measurements in heavy-ion collisions 
at Brookhaven National Laboratory's Relativistic Heavy- 
Ion Collider (RHIC) fll-tlBl]. Both long-range correla- 
tions in pseudo-rapidity and double-peak structures on 
the away-side in A?7 p -A</)-corrclations have been repro- 
duced using initial states with fluctuations in the trans- 
verse plan e, extending as flux-tubes along the beam-line 
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In this work we report on results for both elliptic and 
triangular flow obtained with an event- by- event 3+1D 
relativistic hydrodynamic simulation, an extension of 
MUSIC [7JI , including shear viscosity. We first briefly de- 
scribe the inclusion of viscosity and leave a more detailed 
description to a forthcoming work. 

In the first order formalism for viscous hydrodynamics, 
the stress-energy tensor is decomposed into 



where 



= (e + V)u»u v - Vg^ 



(1) 



(2) 



is the ideal fluid part with flow velocity local energy 
density e and local pressure V . The flow velocity is de- 
fined as the time-like eigenvector of 



Tr v u 



(3) 



with the normalization u v u v = 1 and the pressure is 
determined by the equation of state as a function of e. 

The viscous part of the stress energy tensor in the first- 
order approach is given by 



V 



(4) 



where A^ v = g^ v — u^u v is the local 3-metric and 
V M = A^ iV d v is the local space derivative. Note that 
S^ v is transverse with respect to the flow velocity since 
A^Ui/ = and u v u v = 1. Hence, is also an eigen- 
vector of the whole stress-energy tensor with the same 
eigenvalue e. 



This form of viscous hydrodynamics is conceptually 
simple. However, this Navier-Stokes form is known to in- 
troduce unphysical superluminal si gna ls. There are sev- 



eral remedies for this problem [17H21| . all of them em- 



ploying the second order formalism. In this work, we use 
a variant of the Israel-Stewart formalism derived in [22| , 
where the stress-energy tensor is decomposed as 



j-HV _ rpflV + W ^ _ 

The evolution equations arc 9mT mi/ = and 



(5) 



A£A" g u a d a W af} = -— (W^ - SH - \w^{d a u a ) . 

T n 3 

(6) 

In the r, rj s coordinate system we use, these equations 
can be re-written as hyperbolic equations with sources 



fa J id 



-d a W ab + F b 



and 



d a (u a W cd ) = -{\/T v ){W cd - S cd ) + G 



.:</ 



(7) 



(8) 



where F b and G cd contain terms introduced by the coor- 
dinate change from t, z to r, r\ a as well as those introduced 
by the projections in Eq.(|6|). 

Our approach to solve these hyperbolic equations re- 
lies on the Kurganov-Tadmor (KT) scheme [23j, |24| , to- 
gether with Heun's method to solve resulting ordinary 
differential equations. For details, see [lj|. The main 
difference between the method employed here and the 
method used to solve ideal hydrodynamics in [l6| is 
the appearance of time-derivatives in the source term. 
These are handled with the first order approximation 
g( T n) = (s( T n) — 5 , ( T n-i))/AT in the first step of the 
Heun method, and in the second step we use g{T n ) = 
(9*i T n+i) — g{ T n))l 'At where g*{r n+ i) is the result from 
the first step. 

As in most Eulerian algorithms, ours also suffers from 
numerical instability when the density becomes small 
while the flow velocity becomes large. Fortunately this 
happens late in the evolution of the system. Regularizing 
such instability has no strong effects on the observables 
we are interested in. Some ways of handling this are 
known (for instance see Ref. 25j] ) . 
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FIG. 1: (Color online) Energy density distribution in the transverse plane for one event with 6 = 2.4 fm at the initial time 
(left), and after r = 6 fm/c for the ideal case (middle) and with rj/s = 0.16 (right). 



In this study, we found that setting the local viscosity 
to zero when finite viscosity causes negative pressure in 
the cell as advocated in 25| and reducing the ideal part 
by 5 % works well to stabilize the calculations without 
introducing spurious effects. 

While in standard hydrodynamic simulations with av- 
eraged initial conditions all odd flow coefficients vanish 
by definition, fluctuations generate triangular flow V3 as 
a response to the finite initial triangularity. 

We follow [l5[ and define an event plane through the 
angle 



ib - 1 arctan (Pr sin (™ft)) 
tyn — arcran - - — — , 

n (j>t cos(n(p)) 



(9) 



where the weight px is chosen for best accuracy [2f 
Then, the flow coefficients can be computed using 



v n = (cos(n((/>- ip n ))) . 



(10) 



The initialization of the energy density is done using 
a Glauber Monte-Carlo model (see Before the col- 

lision the density distribution of the two nuclei is de- 
scribed by a Woods-Saxon parametrization, which we 
sample to determine the positions of individual nucleons. 
The impact parameter is sampled from the distribution 
P(b)db = 2bdb/{b 2 aaK -b 2 a{n ), where & min and 6 max depend 
on the given centrality class. Then we determine the dis- 
tribution of binary collisions and wounded nucleons. Two 
nucleons are assumed to collide if their relative transverse 
distance is less than D = a nn/tt, where <jnn is the in- 
elastic nucleon-nucleon cross-section, which at top RHIC 
energy of y/s = 200AGcV is (Jnn — 42 mb. The energy 
density is distributed proportionally to the wounded nu- 
cleon distribution. For every wounded nuclcon we add a 
contribution to the energy density with Gaussian shape 
(in x and y) and width ao — 0.4 fm. In the rapidity 
direction, we assume the energy density to be constant 
on a central plateau and fall like half-Gaussians at large 
I % I (see [HI). This procedure generates flux-tube like 
structures compatible with measured long-range rapidity 
correlations 28 - 3(J • The absolute normalization is deter- 
mined by demanding that the obtained total multiplicity 
distribution reproduces the experimental data. 



As equation of state we employ the parametrization 
"s95p-vl" from (3l| . obtained from interpolating between 
lattice data and a hadron resonance gas. 

In Fig.[TJ we show the energy density distribution in 
the transverse plane for an event with impact parameter 
b = 2.4 fm at the initial time r = 0.4 fm/c and at time 
r = 6 fm/c for rj/s = and rj/s = 0.16. This clearly 
shows the effect of dissipation. 

We perform a Cooper-Frye freeze-out using 



E 



dN 



dN 



d 3 p dypTdpTC 



r = 9i , (11) 

5 p Js 



where gi is the degeneracy of particle species i, and E 
the freeze-out hyper-surface. In the ideal case the distri- 
bution function is given by 



/Kjv) = /o(u^) 



1 



1 



(2tt) 3 exp(( U v Pll - im)/Tfo) ± 1 



(12) 

where [ii is the chemical potential for particle species 
i and Tfo is the freeze-out temperature. In the finite 
viscosity case we include viscous corrections to the dis- 
tribution function, / = /o + 5f, with 



Sf = f (l ± f )p a pPW t 



1 



aft 



2(e + T)T 2 



(13) 



where W is the viscous correction introduced in Eq. (|5|) . 
Note that the choice Sf ~ p 2 is not unique [HJ . 

The algorithm used to determine the freeze-out surface 
E has been presented in [l6|. It is very efficient in de- 
termining the freeze-out surface of a system with fluctu- 
ating initial conditions. To demonstrate this, we present 
the freeze-out surface in the a;-r-plane in the vicinity of 
y = fm and r] s = for two different initial distribu- 
tions compared to that for an averaged initial condition 
in Fig. [2] The arrows are projections of the normal vector 
on the hyper-surface element onto the x-t plane. 

We include resonances up to the </>-meson. We found 
that the pseudorapidity dependence of both V2 and V3 is 
affected notably by the inclusion of resonance decays, im- 
proving the agreement of l>2(?7p) with data significantly. 
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2: (Color online) Freeze-out surfaces for two different 




events (red and yellow) compared to that for the averaged 
initial condition (gray). 
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FIG. 3: (Color online) Charged hadron «2 for different cen- 
tralizes as a function of transverse momentum for averaged 
initial conditions (avg) and event-by-event simulations (e-b-e) 
using different viscosity to entropy density ratios compared to 
STAR [a3 and PHENIX [H data. 

V2(pt) a t midrapidity is almost unaffected by the reso- 
nances while Vs(pt) is reduced by approximately 20-30%. 

Fig.[3]shows the elliptic flow V2 for charged hadrons as a 
function of transverse momentum obtained from an aver- 
aged initial condition in the ideal case and for an average 
over 100 individual events forn/s e {0,0.08,0.16}. We 



341. The 



compare to data from STAR |33j and PHENIX 
used minimal, maximal, and average impact parameter 
for each centrality class are given in Tabic fl] . While in 
the most central collisions fluctuations increase V2 com- 



centrality [%] 


femin [fm] 


&max [fm] 


(b) [fm] 


0-5 





3.37 


2.24 


10-20 


4.75 


6.73 


5.78 


15-25 


5.83 


7.53 


6.7 


30-40 


8.23 


9.5 


8.87 



TABLE I: Used impact parameters. 

pared to the case with averaged initial conditions, for 10- 
20% central collisions the difference is negligible and for 
30-40% central collisions fluctuations reduce the elliptic 
flow. The increase for central collisions is easy to un- 
derstand since we are now determining i>2 in every event. 
Single events have a larger anisotropy with respect to the 
event-plane than the average with respect to the reaction 
plane, hence increasing the obtained V2- This effect de- 
creases with increasing centrality eventually making the 
event-by-event i>2 smaller compared to the averaged ini- 
tial condition case. This can be understood by the fact 
that for more peripheral collisions, lumps in the initial 
condition tend not to align perfectly with the statisti- 
cally determined event plane. 

Viscosity reduces the elliptic flow for all centralities as 
also found in (2+l)-dimensional simulations [35 - 38 1 . 

Triangular flow V3 as a function of transverse momen- 
tum is shown in Fig.|4] V3 depends less strongly on 
the centrality than V2 since it is completely fluctuation 
driven. It is largest for an ideal fluid and reduces simi- 
larly to v 2 with increasing viscosity of the medium. 

The upper panel of Fig.[5]shows the pseudo-rapidity de- 
pendence of V2 for 15-25% central collisions compared to 
PHOBOS data (39|. A reduction of elliptic flow with in- 
creasing viscosity is visible, particularly for large pseudo- 
rapidities |?7p|, which has been anticipated Til 3. In the 
lower panel of Fig. [5] we present the ^-dependence of 1)3. 
Again, the decrease of V3 with increasing viscosity is vis- 
ible, being strongest for large \rj p \. 

We presented the elliptic and triangular flow coef- 
ficients obtained with an event-by-event analysis us- 
ing (3+l)-dimcnsional relativistic viscous hydrodynam- 
ics. Charged hadron elliptic flow around midrapidity is 
well described for a wide range of centralities when using 
77/s = 0.08, the conjectured lower bound from AdS/CFT 
41| . A similarly small value was found in a parton cas- 
cade model based on perturbative QCD [HI]. Larger vis- 
cosities underestimate elliptic flow. Shear viscosity re- 
duces V2 especially for larger pseudo-rapidities, however, 
the data is still overestimated away from midrapidity. 
Triangular flow has a weaker dependence on centrality. 
We determined its transverse momentum and pseudo- 
rapidity dependence, as well as its dependence on n/s. 
When triangular flow data becomes available, combined 
analyses of both V2 and V3 can make an accurate deter- 
mination of the shear-viscosity possible. 
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FIG. 4: (Color online) Charged hadron vs for different cen- 
tralities as a function of pr for event-by-event simulations 
using different viscosity to entropy density ratios. 
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FIG. 5: (Color online) Charged hadron vi and vz as a func- 
tion of pseudo-rapidity for event-by-event simulations using 
different viscosity to entropy density ratios. 
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